Implementation of the Conjugate GradientAlgorithm

Computation of the inner state parameters in DSO inversion requires solving a large normal matrix system. A combined conjugate gradient and Lanczos iterative technique can be used to both solve the system and approximate some of the spectrum of the normal operator. At each iteration of the conjugate gradient algorithm, a small tridiagonal matrix (of dimension equal to the number of iterations) is created which has extreme eigenvalues approximating those of the original matrix. A second matrix whose columns are the normalized residual vectors from the conjugate gradient algorithm allows the corresponding eigenvectors to be computed as well if desired. Implemented so that it can be applied to diierent DSO inversion problems, the conjugate gradient code provides the user with a tool to analyze the condition of the problem as well as the quality of the inversion results. Storage of the residual (Lanczos) vectors may be costly if large problems are being solved. Numerical experiments indicate that the largest eigenvalue is the best approximation in the spectrum and the smallest is generally the next best.